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Abstract 


Classically,  fluid  dynamics  have  been  dealt  with  analytically  because  of  the  lack  of  nu¬ 
merical  resources  (Yeung,  1982).  With  the  development  of  computational  ability,  many 
formulations  have  been  developed  which  typically  use  the  traditional  Navier-Stokes  equa¬ 
tions  along  with  an  Eulerian  grid.  Today,  there  exists  the  possibility  of  using  a  moving  grid 
(Lagrangian)  along  with  a  meshless  discretization. 

The  first  issue  in  meshless  fluid  dynamics  is  the  equations  of  motion.  There  are  currently 
two  types  of  Lagrangian  formulations.  Spherical  Particle  Hydrodynamics  (SPH)  is  a  method 
which  calculates  all  equations  of  motion  explicitly.  The  Moving  Particle  Semi-implicit  (MPS) 
method  uses  a  mathematical  foundation  based  on  SPH.  However,  instead  of  calculating 
all  laws  of  motion  explicitly,  a  fractional  time  step  is  performed  to  calculate  pressure.  A 
proposed  method,  Lagrange  Implicit  Fraction  Step  (LIFS),  has  been  created  which  improves 
the  mathematical  formulations  on  the  fluid  domain.  The  LIFS  method  returns  to  Continuum 
mechanics  to  construct  the  laws  of  motion  based  on  decomposing  all  forces  of  a  volume. 
It  is  assumed  that  all  forces  on  this  volume  can  be  linearly  superposed  to  calculate  the 
accelerations  of  each  mass.  The  LIFS  method  calculates  pressure  from  a  boundary  value 
problem  with  the  inclusion  of  proper  flux  boundary  conditions. 

The  second  issue  in  meshless  Lagrangian  dynamics  is  the  calculation  of  derivatives  across 
a  domain.  The  Monte  Carlo  Integration  (MCI)  method  uses  weighted  averages  to  calculate 
operators.  However,  the  MCI  method  can  be  very  inaccurate,  and  is  not  suitable  for  sparse 
grids.  The  Radial  Basis  Function  (RBF)  method  is  introduced  and  studied  as  a  possibil¬ 
ity  to  calculate  meshless  operators.  The  RBF  method  involves  a  solution  of  a  system  of 
equations  to  calculate  interpolants.  Machine  expenses  are  shown  to  limit  the  viability  of 
the  RBF  method  for  large  domains.  A  new  method  of  calculation  has  been  created  called 
Multi-dimensional  Lagrange  Interpolating  Polynomials  (MLIP).  While  Lagrange  Interpolat¬ 
ing  Polynomials  are  essentially  a  one-dimensional  interpolation,  the  use  of  “dimensional-cuts” 
and  Gaussian  quadratures  can  provide  multi-dimensional  interpolation. 

This  paper  is  divided  into  three  sections.  The  first  section  specifies  the  equations  of 
motion.  The  second  section  provides  the  mathematical  basis  for  meshless  calculations.  The 
third  section  evaluates  the  effectiveness  of  the  meshless  calculations  and  compares  two  fluid- 
dynamic  codes. 
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Chapter  1 

Equations  of  Motion 


Eulerian  fluid  dynamics  calculates  the  flow  of  fluid  through  a  given  volume.  Lagrangian 
dynamics  calculates  the  movement  of  discrete  volumes  in  a  flow.  We  will  utilize  Lagrangian 
forms  of  equations  to  model  highly  deformable  flows,  while  tracking  the  movement  of  discrete 
volumes.  All  types  of  Lagrangian  formulations  will  use  a  variation  of  the  Navier-Stokes 
equations. 


1.1  Navier-Stokes  Equations 


Continuum  mechanics  provides  a  context  to  construct  the  laws  of  motion.  We  will  use 
balance  laws  of  linear  momentum,  angular  momentum,  and  mass  (Chadwick,  1999).  The 
basic  postulate  for  linear  momentum  on  a  material  volume  is  from  Newton’s  second  law 
(Johnson,  1990): 


There  exists  a  frame  of  reference  for  which,  at  any  instant,  the  rate  of  change  of 
linear  momentum  of  a  material  volume,  fl,  is  equal  to  the  resultant  force  acting 
on  the  mass  in  that  volume. 


=  =  U!) 

in  which  L  is  a  force  vector,  dQ  is  the  surface  of  the  material  volume,  ||  is  the  material 
derivative,  p  is  the  density  of  the  fluid,  and  u  is  the  velocity  vector. 

Cauchy  developed  the  idea  of  a  stress  tensor,  u,  in  which: 


Can  — 


a  ■  ndS 


Cn  = 


an 


fdn 


(1.2) 


where  n  is  the  outward  unit  normal  and  /  are  the  body  forces.  By  use  of  the  divergence 
theorem  we  can  relate  the  surface  integral  to  a  volume  integral. 


a ■ ndS  = 


an 


V-cxdfi 


(1.3) 
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Therefore,  our  balance  of  linear  momentum  can  be  expressed  as: 

MIL  «d!i  =  III  v'gdtt+ III /dn  (L4) 

Dealing  with  only  volume  integrals,  the  integrand  must  vanish  within  the  fluid,  and  we  get 
Cauchy’s  equation  of  motion: 


^-pu  =  V  ■  g  +  l  (1.5) 

To  locally  conserve  angular  momentum,  we  specify  that  the  stress  tensor  is  symmetric. 

a  =  aT  (1.6) 

Finally,  we  will  specify  conservation  of  mass. 

<i7) 

Using  the  transport  theorem  and  specifying  that  the  integrand  is  zero  throughout  the  fluid, 
we  can  express  conservation  of  mass  as: 

— ^  +  pV  •  u  =  0  (1.8) 


1.1.1  Constitutive  relations 


We  will  specify  a  Newtonian  constitutive  relation  (Newman,  1977). 


<7=  —pl  +  fiD  (1.9) 

in  which  p  is  the  pressure,  /  is  the  identity  tensor,  //  is  the  viscosity  coefficient,  and  D  is  the 
symmetric  deformation  tensor. 


-p 

0 

0 

9  du 
^  dx 

du  |  dv 
dy  '  dx 

du  _i_  dw  - 
dz  '  dx 

>]  = 

0 

-p 

0 

+  /i 

dv  _i_  du 
dx  dy 

9  dv 
dy 

dv  _i_  dw 
dz  '  dy 

0 

0 

- p  _ 

dw  |  du 
-  dx  '  dz 

dw  |  dv 
dy  '  dz 

9  dw 

Z  dz  J 

(1.10) 


1.1.2  Boundary  Conditions 

Our  fluid  domain  (D)  is  shown  in  figure  ED-  There  are  two  boundary  conditions  on  the 
domain.  Dirichlet  boundary  conditions  apply  to  the  free  surface.  Neumann  boundary  con¬ 
ditions  apply  to  the  surface  next  to  the  wall.  The  walls  are  not  a  part  of  the  fluid  domain, 
but  serve  to  administer  the  Neumann  boundary  conditions. 
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A  diagram  of  the  domain 

•  The  fluid  domain  fi  is  shown  in  blue 

•  There  are  Neumann  boundaries  Tt  and 
Dirichlet  boundaries 

•  There  is  closure  such  that  =  Tt  U  T^ 

•  Outward  unit  normals  are  shown  as  n 


Figure  1.1:  Domain  used  in  equations  of  motion 


1.1.3  Lagrangian  Forms 

Part  of  the  ease  of  Lagrangian  coordinates  is  being  able  to  keep  the  material  derivatives 
instead  of  using  a  convection  term.  In  doing  so,  we  express  the  equations  of  motion  for 
discrete  volumes  of  fluid.  These  individual  volumes  move  with  the  flow  of  the  fluid.  If  we 


1995]) 


enforce  water  fluid  properties  such  that  ^  =  0,  then  the  law  of  motion  is  simply:  (Brodkey, 


where  a  is  the  acceleration  vector  which  is  equal  to  the  material  derivative  of  velocity.  If  we 
use  explicit  time-stepping,  we  can  obtain: 


Su  =  —  (V  •  a  +  /) 
P  ~ 

in  which, 


un+1  =  u11  +  5u 


(1.12) 


(1.13) 


Position  vectors  are  therefore: 


xn+1  =  xn  +  5tun+l 


The  continuity  law  is: 


—  P  (V  •  u)  =  0 


(1.14) 


(1.15) 
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1.2  Spherical  Particle  Hydrodynamics 

Monaghan  first  suggested  the  use  of  individual  particles  in  the  calculation  of  fluid  dynamics. 


Each  particle  interacts  with  neighboring  particles  according  to  the  law  of  motion.  (Mon¬ 


aghan,  1994) 


SPH  is  calculated  in  a  fully  explicit  manner.  Using  eqn(l.ll),  we  calculate  a  explicitly 
as  a  function  of  particles  a  and  b: 


_  Pa  Pb  yr 

~  ^ab  ~  +  3  + 

Pa  Pb 


(1.16) 


Pressure  is  calculated  as  an  equation  of  state  for  a  very  “stiff”  gas.  It  is  a  function  of  the 
density  calculated  at  a  time  step  n  versus  the  standard  density  (p°): 


(1.17) 


The  parameter  7  is  set  arbitrarily  at  7,  and  the  coefficient  B  is  somewhat  defined  by  the 
speed  of  sound.  Density  is  calculated  for  each  particle  as: 


Pa=  madtta 


(1.18) 


Both  bulk  and  shear  viscosity  are  calculated  through: 


n 


ab 


-acpab  +  /fyqfe 

1/2  {pa  +  Pb) 


(1.19) 


The  parameter  fj,ab  roughly  tries  to  simulate  the  laplacian  operator,  a  introduces  the  viscosity, 
c  is  the  speed  of  sound,  and  (3  is  usually  set  to  zero. 

Boundary  conditions  for  SPH  are  treated  in  an  ad-hoc  approach.  Boundary  particles  are 
treated  no  differently  than  any  other  particle,  regardless  of  if  it  is  on  the  free  surface  or  is 
next  to  a  wall.  The  most  common  way  to  ensure  that  particles  in  SPH  do  not  violate  wall 
boundaries  is  through  an  extra  force.  When  a  particle  comes  near  a  wall,  a  repulsive  force 
pushes  on  the  particle  such  that  it  does  come  into  the  wall.  The  degree  of  force  required, 
however,  must  be  set  for  different  types  of  simulations. 

The  basic  algorithm  for  this  method  is  shown  in  table  |1.1[ 

There  have  been  very  good  results  with  the  SPH  equations  including  dam  collapse  and 
breaking  waves  (Monaghan,  1994).  Furthermore,  SPH  has  been  shown  to  be  robust  in 
“floating  object”  problems  including  a  falling  wedge  penetrating  water  as  well  as  a  wave 


interaction  with  a  floating  box  (Doring  et  ah,  2002). 
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INITIALIZE  the  domain 

CALCULATE  the  initial  density,  (p°) 


eqn(1.18) 


FOR  EACH  time  step,  n 

CALCULATE  the  density,  ( pn ) 
EVALUATE  the  pressure,  (pn) 
EVALUATE  the  viscosity,  (IT1) 
EVALUATE  the  motion,  8u 
UPDATE  the  velocity,  ( un+1 ) 
UPDATE  the  position,  ( xn+1 ) 

END 


eqn' 

eqn' 

eqn' 

eqn' 

eqn' 

eqn' 


(1-18) 


(1-17) 


(1-19) 


(1.12) 


(1.13) 


(1.14) 


Table  1.1:  SPH  Algorithm 


1.3  Moving  Particle  Semi-Implicit 


The  MPS  equations  are  a  variation  on  SPH  in  that  instead  of  pressure  being  calculated 
explicitly,  it  is  solved  via  Poisson’s  equation  (Yoon  et  ah,  1999).  This  is  done  by  a  fractional 
step  which  splits  the  laws  of  motion  into  two  parts  and  uses  the  continuity  equation  to 
connect  the  two. 

The  first  fractional  step  is  defined  as: 


u*  =  un  +  —  (vX72un  +  pg) 

p 


(1.20) 


The  second  fractional  step  is  defined  as: 


u**  =  Vpn+1 
P 


(1.21) 


In  order  to  calculate  the  second  fractional  step,  we  must  find  a  way  to  calculate  the 
pressure.  This  can  be  done  by  first  finding  the  change  in  density  due  to  the  first  fractional 
step 


dp*  _  1  p*  -  p° 

dt  St  p° 


(1.22) 


Pressure  is  calculated  as  a  solution  to  Poisson’s  equation.  There  is  only  one  boundary 
condition  in  which  particles  on  the  free  surface  have  Dirichlet  conditions  where  the  pressure 
is  set  to  zero.  Particles  near  walls  do  not  have  a  boundary  condition,  but  are  instead  treated 
by  the  held  equation.  The  right  hand  side  of  the  held  equation  is  the  divergence  of  the  hrst 
fractional  velocity: 


vV+1  = 


P_[ dp 

St  dt 


(1.23) 
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Finally,  the  total  velocity  is  the  summation  of  fractional  velocities: 


un+l  =  u* 


i  * 

U 


(124) 


Boundary  conditions  in  MPS  have  more  treatment  than  in  SPH.  First,  particles  on  the  free 
surface  have  the  pressure  set  to  zero.  Particles  next  to  the  walls,  however,  do  not  have  any 
boundary  conditions.  Instead,  the  walls  are  treated  as  fluid  particles  which  are  stationary. 
This  procedure  is  easy  to  implement,  however  it  violates  the  equations  of  motion.  The  laws  of 
motion  are  expressed  in  terms  of  integrals  over  a  fluid  domain  as  well  as  material  derivatives. 
When  a  derivative  is  calculated  which  is  not  material,  than  the  equation  of  motion  is  invalid. 


The  basic  algorithm  for  this  method  is  shown  in  table  1.2 


INITIALIZE  the  domain 

CALCULATE  the  initial  density,  (p°) 


eqn(1.18) 


FOR  EACH  time  step,  n 

EVALUATE  the  first  fractional  velocity,  («*) 
CALCULATE  the  density,  (p*) 

CALCULATE  the  change  in  density,  (^§f) 

SOLVE  for  pressure,  ( pn+1 ) 

EVALUATE  the  second  fractional  velocity,  (u**) 
SUMMATE  the  fractional  velocities,  (un+1) 
UPDATE  the  position,  (:rn+1) 

END 


(1. 


eqm 
eqni 
eqni 
eqni 
eqn(l 


(1.20) 


(1 


eqn(l 

eqn(l 


(1 


18) 

22) 

23) 

2T) 

24) 
14) 


Table  1.2:  MPS  Algorithm 


Again,  good  results  have  been  obtained  for  a  sloshing  tank  (Yoon  et  ah,  1999),  breaking 
waves  (Lo  and  Shao,  2002),  and  wave  generation  (Gotoh  and  Sakai,  1999).  However,  it 
should  be  noted  that  the  differential  equation  (1.23|)  only  has  boundary  conditions  on  one 
side  of  the  fluid  and  is  therefore  an  ill-posed  boundary  value  problem. 


1.4  Lagrange  Implicit  Fraction  Step 

The  LIFS  equations  were  created  to  use  the  principles  of  SPH  and  MPS,  but  to  improve  the 
mathematical  foundations.  For  each  force  term  in  the  balance  of  linear  momentum,  we  break 
it  into  a  corresponding  explicit  and  implicit  velocities.  We  then  solve  the  continuity  equation 
for  the  implicit  velocities  which  ensures  that  the  momentum  equation  does  not  result  in  a 
change  in  density. 

To  reiterate  our  explicit  procedure,  we  take  the  Lagrangian  conservation  of  momentum 
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and  break  it  into  discrete  time  steps: 


~P  (V  -  2  +  /) 

(1.25) 

—  (V  '  (7  +  /) 

P  - 

(1.26) 

un  +  Su 

(1.27) 

un+1  = 

We  assume  that  the  constitutive  relation  eqn(|1.9[)  and  our  summation  of  forces  is  linear: 

a  —  -pL+  pH  (1-28) 

a  a 

For  each  time  step,  we  will  define  intermediate  velocities  based  on  the  forces  which  con¬ 
tributed  to  them: 


Su  =  Su  +  Su  +  Su 


(1.29) 


The  tilde  velocity  is  the  motion  due  to  body  forces,  the  hat  velocity  is  the  motion  due  to 
viscosity,  and  the  spherical  velocity  is  the  motion  due  to  pressure. 


Su  = 

(1.30) 

P~ 

Su  = 

St,  A 

— V  •  a 

(1.31) 

P 

Su  = 

St,  0 
—  V  •  g_ 

(1.32) 

P 

In  each  time  step,  we  have  both  implicit  and  explicit  parts.  We  define  the  body  forces 
explicitly  and  due  only  to  gravity: 

Su  =  (St)  (g)  (1.33) 

The  hat  velocity  is  explicitly  defined  from  viscosity  calculated  from  the  current  velocities: 

Su  =  (St)  i/VV  (1.34) 


The  spherical  velocity  is  calculated  implicitly  from  the  pressure  (Yeung  and  Ananthakrish- 


nan,  1992): 


Su  =  — — Vpn+1 
P 

We  can  obtain  the  pressure  through  the  continuity  equation, 


(1.35) 


^  +  pV  •  un+1  =  0 
Vt  H  ~ 


Vp 


+  pV  •  (un  +  Su  +  Su  +  Su)  =  0 
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(1.36) 

(1.37) 


Assuming  p  is  constant  in  time,  we  separate  the  explicit  and  implicit  velocity  terms  (Anan- 
thakrishnan  and  Yeung,  1994): 


V  •  Su  =  —V  •  ( un  +  8u  +  Su )  (1.38) 

By  applying  our  constitutive  relation  into  our  spherical  velocity,  the  continuity  equation 
becomes: 

V  '  (7  (V^n+1))  =  V  •  (un  +  Su  +  Su)  (1.39) 

Because  we  are  working  with  incompressible  fluid,  we  can  show  that: 

V-un  =  0  (1.40) 


Furthermore,  since  the  body  forces  are  simply  a  constant  scalar  field: 

V  •  Su  =  0 


(1.41) 


Therefore,  the  only  term  left  is  due  to  viscosity.  This  can  be  found  to  be  Poisson’s  equation: 

**"-&■*-*%  (1'42) 

This  can  be  setup  as  a  boundary  value  problem  by  making  specifications  on  the  field  equations 
and  the  boundary  conditions: 


(1.43) 


The  third  prescription  can  be  found  through: 


1)  V2pn+1  =  t-V  •  Su  in  n 

Ot 

2)  pn+1  =  0  on  Td,  the  free  surface 

3)  un+1  ■  n  =  0  on  Tt,  the  walls 


un+1 • n 


(; un  +  8u  +  8u  +  Su)  ■  n 
Su  ■  n 


0 

0 

—  ( un  +  Su,  +  Su)  •  n 


(1.44) 


Now,  we  insert  eqn(1.35)  into  eqn(1.44). 
direction  as:  Vn. 


We  express  the  directional  derivative  in  the  n 


Vnpn+1  =  (un  •  n  +  Su  ■  n  +  Su  ■  n) 
St 


(1.45) 


We  can  also  simplify  this  equation  in  that  yJ1  •  n  =  0.  Therefore,  our  boundary  condition  is: 


VnA 


,n+ 1  _  P_ 

St 


(Su  ■  n  +  Su  ■  n) 


(1.46) 


An  algorithm  for  the  Lagrange  Implicit  Fraction  Step  method  is  shown  in  table  1.3 


INITIALIZE  the  domain 

FOR  EACH  time  step,  n 

EVALUATE  the  body  forces  for  Su  eqn1 

EVALUATE  the  viscous  forces  for  Su  eqn1 
EVALUATE  the  divergence  for  the  BVP  eqn' 
EVALUATE  the  flux  BCs  for  the  BVP  eqn' 
SOLVE  Poisson’s  equation  for  pn+1  eqn' 

EVALUATE  the  spherical  forces  for  Su  eqn' 
EVALUATE  the  motion,  (Su)  eqn' 

UPDATE  the  velocity,  (un+1)  eqn' 

UPDATE  the  position,  (xn+1)  eqn> 

END 


(1.33) 


(1.34) 


(1.42) 


(1.46) 


(1.43) 


(1.35) 


(1-29) 


(1.13) 


(1.14) 


Table  1.3:  LIFS  Algorithm 

An  interesting  aspect  is  when  dealing  with  an  inviscid  fluid.  When  this  occurs,  Su  =  0. 
Therefore,  the  field  equation  becomes  Laplace’s  equation: 


vV+1  =  o 

Furthermore,  the  flux  boundary  conditions  are  only  a  function  of  body  forces: 


Vnp' 


.71+ 1  P_ 


—5u  •  n 
ot 


(1.47) 


(1.48) 


If  we  set  viscosity  to  zero,  we  can  simplify  the  algorithm,  shown  in  table  1.4 


INITIALIZE  the  domain 

FOR  EACH  time  step,  n 

EVALUATE  the  body  forces  for  Su  eqn' 

EVALUATE  the  flux  BCs  for  the  BVP  eqn' 

SOLVE  Laplace’s  equation  for  pn+l  eqn> 

EVALUATE  the  spherical  forces  for  Su  eqn' 
EVALUATE  the  motion,  (Su)  eqn' 

UPDATE  the  velocity,  («n+1)  eqn' 

UPDATE  the  position,  (xn+1)  eqn' 

END 


(1.33) 


(1.48) 


(1.47) 


(1.35) 


(1.29) 


(1-13) 


(1-14) 


Table  1.4:  LIFS  Inviscid  Algorithm 
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Chapter  2 

Methods  of  Computation 


After  the  equations  of  motion  are  set,  the  problem  becomes  how  to  calculate  the  equations 
between  particles.  Because  all  the  methods  are  “meshless,”  problems  of  constructing  grids 
are  alleviated.  However,  the  question  of  how  particles  interact  still  remain. 

The  laws  of  motion,  for  a  meshless  method,  are  extremely  dependent  on  how  the  par¬ 
ticles  interact.  In  essence,  the  methods  of  computation  represent  an  interpolation  of  the 
solution  between  particles.  The  better  the  interpolation,  the  better  the  solution  matches  the 
appropriate  laws. 

For  every  method  of  computation,  what  we  wish  to  find  is  a  continuous  function  f(x)  as 
a  function  of  discrete  values,  f(x  =  xt)  =  fix,)-  Methods  of  computation  typically  have  the 
form  of  finding  a  continuous  function  by  multiplying  a  “kernel”  (w)  by  a  value  (a)  evaluated 
at  each  discrete  location. 


n—  1 

f(x)  =  J2wl(rJ)a(xJ)  (2.1) 

j=o 

These  methods  use  a  nodal-centered  position  vector  as  an  important  value  in  the  weighting 
function.  We  will  define  the  radius,  ry ,  as  the  Euclidean  norm  of  two  position  vectors: 

rv  =  (2-2) 

Most  particle  methods  evaluate  a  function  g  with  the  origin  centered  at  the  particle  location 
xt.  The  argument  for  the  function  is  the  distance  from  another  point  £  to  the  particle  center. 
We  will  define  the  notation  for  this  function  g  as  a  function  of  the  radius: 


9^i)  =  9{rtJ)=gt{rj)  (2.3) 

All  of  these  methods  of  computation  are  discretely  calculated  via  the  use  of  matrices.  For 
each  method  of  computation,  the  desired  properties  of  each  method  includes  at  least  the 
formation  of  a  matrix  which  can  perform  interpolation,  (w),  as  well  as  matrices  which  can 
calculate  gradients  ( dwx  and  dwy ).  Furthermore,  both  the  MPS  and  the  LIFS  equations 
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require  the  calculation  of  the  Laplacian  operator.  This  operator  can  be  represented  in  matrix 
form  as  ddw. 

Every  one  of  these  matrices  required  by  the  equations  of  motion  will  be  constructed 
through  the  use  of  the  following  methods.  The  Monte  Carlo  Integration  (MCI)  method 
calculates  both  the  interpolation  and  derivatives  based  on  a  process  of  averaging  neighbors. 
The  Radial  Basis  Function  (RBF)  method  calculates  these  operators  by  first  solving  a  linear 
system  of  equations.  After  one  matrix  is  inverted,  all  differencing  operations  can  be  per¬ 
formed.  Multi-dimensional  Lagrange  Interpolating  Polynomials  (MLIP)  calculates  operators 
in  a  meshless  manner  by  creating  Lagrange  Interpolating  Polynomials  and  their  derivatives. 


2.1  Monte  Carlo  Integration 

The  Monte  Carlo  Integration  method  is  the  traditional  method  for  calculating  interpolants 
in  a  meshless  manner.  The  principal  behind  the  MCI  method  is  that  information  about  one 
particle  is  stored  in  its  neighbors.  Therefore,  we  can  use  a  weighted  average  of  the  neighbor 
particles  to  not  only  obtain  the  interplant,  but  also  to  find  derivatives. 

2.1.1  Formulation 

By  using  an  interpolating  weighting  function,  discrete  values  can  be  found  from  the  neighbors 
in  the  fluid  volume  (f2).  The  weighting  kernel  w  at  each  node  i  is  a  function  of  the  distance, 
r3,  and  the  extent  of  the  kernel,  h.  The  extent  of  the  kernel,  h,  determines  how  many  other 
particles  are  used  in  the  calculations.  We  will  define  h  as  a  resolution  parameter.  As  the 
resolution  parameter  increases,  more  neighbor  particles  are  used  in  the  calculations,  and  the 
stability  of  the  calculations  increases.  As  the  resolution  parameter  decreases,  less  neighbor 
particles  are  used,  and  the  accuracy  of  the  calculations  increases. 

f(x)=  [  f(§)wt(r3,h) dfi  (2.4) 

Jn 

2.1.2  Prescriptions  on  the  kernel 

There  are  three  prescriptions  on  the  weighting  kernel.  The  first  is  that  the  function,  w,  is 
normalized,  such  that  energy  is  not  added  nor  subtracted  from  the  system. 

I  wl{r3 ,  h)d£l  =  1  (2.5) 

Ja 

The  second  prescription  is  that  the  weighting  function  must  become  the  Dirac  delta  function 
as  the  interaction  radius,  h,  becomes  small: 


A  diagram  of  Monte  Carlo  Integration 

•  The  red  particle  is  the  particle  to  which 
integration  is  being  performed 

•  The  white  particles  are  neighbor  parti¬ 
cles 

•  The  length  h  represents  the  maximum 
interaction  of  the  center  particle  on  the 
surrounding  particles 


Figure  2.1:  Monte  Carlo  Integration 

The  third  prescription  is  that  the  weighting  function  has  a  finite  radius: 

lim  wt(\rj\,h)  =  0  (2.7) 

b \^h 

This  third  prescription  allows  us  to  redefine  our  limits  of  integration.  Instead  of  the  entire 
fluid  domain  (fi),  we  can  redefine  the  integration  to  take  place  only  in  the  domain  of  each 
“particle”,  (f2p).  Furthermore,  if  our  kernel  is  smooth  and  continuous  we  can  take  the 
derivative  of  /.  Because  of  finiteness,  we  find  that  the  operation  is  passed  only  to  the  kernel 
since  w(h,  h )  =  0. 

Vf{x)=  f  mVwt{r3,h) dnp  (2.8) 


2.1.3  Evaluation 

Weighting  kernels  typically  have  compact  support  and  are  smooth  to  various  degrees  of 
differentiability.  Therefore,  kernels  are  mostly  based  on  either  splines  or  Gaussian  curves. 
We  will  choose  a  variation  on  a  Gaussian  kernel  in  which: 

w,(r’h)  =  ^exp  (2/?)  +  2  0  +  erf  (yi*))  <2'9) 

V.«.(r,  *)  =  £(: l+«rf(^)) 

1  /— r2\ 

=  wTxpljW 
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(2.10) 

(2.11) 


V2w(r,  h) 


We  can  implement  the  MCI  method  using  three  steps.  The  first  step  is  to  initialize  the 
kernel  matrices.  The  second  step  is  to  use  the  kernel  to  calculate  derivatives.  The  final  step 
is  to  set  up  a  matrix  to  invert  and  solve  a  linear  system. 

Table  12.11  shows  the  initialization  routine  for  the  MCI  method.  All  of  the  matrices  for 
calculating  operators  are  created  here:  interpolation  (W),  x-gradient  (dWx),  y-gradient 
(< dW y ),  and  laplacian  (ddW). 


FOR  i  =  0  to  (n-1) 

FOR  j  =  0  to  (n-1) 
CALCULATE  r[i,j] 
CALCULATE  W[i,j] 
CALCULATE  dWx[i,j] 
CALCULATE  dWy[i,j] 
CALCULATE  ddW[i,j] 
NEXT  j 
NEXT  i 


Table  2.1:  MCI  Initialization  Algorithm 


Table  |2.2|  shows  how  derivatives  of  a  nodally-specified  function,  /,  may  be  calculated 
with  the  operators  (dW x)  and  ( dW y ).  If  the  kernel  was  not  normalized,  we  can  fix  this  by 
dividing  the  resulting  vectors  by  the  sum  of  the  absolute  value  of  the  kernel. 


FOR  i  =  0  to  (n-1) 

FOR  j  =  0  to  (n-1) 

dfdxti]  =  dWx[i,j]  *  f[j]  eqnjB 
dfdy[i]  =  dWy[i,j]  *  f[j]  eqn(2.8) 
NEXT  j 
NEXT  i 


Table  2.2:  MCI  Derivative  Calculation  Algorithm 


Table  |2.3|  gives  an  algorithm  of  how  a  boundary  value  problem  can  be  solved  with  the 
MCI  method.  This  algorithm  is  used  to  solve  equations  such  as  eqn(1.43).  This  linear  system 
of  equations  has  the  form  of  multiplying  a  matrix  ([K])  by  the  solution  vector  ({a})  to  obtain 
the  forcing  vector  ({/}): 


[K]  M  =  {/} 


(2.12) 


Before  calculations  begin,  each  particle  must  be  identified  as  being  a  member  of  the  field 
equation  or  a  boundary  particle.  Boundary  particles  are  then  separated  for  having  Dirichlet 
or  Neumann  boundary  conditions.  Particles  having  Neumann  conditions  have  two  values 
associated  with  them  indicating  the  normals  to  the  wall  in  the  x  and  y  directions  (nx  and 
ny). 
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FOR  i  in  the  field  equation 
FOR  j  =  0  to  (n-1) 

K[i,j]  =  ddW[i,  j] 

NEXT  j 
NEXT  i 

FOR  i  having  Dirichlet  conditions 
FOR  j  =  0  to  (n-1) 

K[i,  j]  =  W[i ,  j] 

NEXT  j 
NEXT  i 

FOR  i  having  Neumann  conditions 
FOR  j  =  0  to  (n-1) 

K[i,j]  =  dWx[i,j]  *  nx[i]  +  dWy[i,j]  *  ny[i] 
NEXT  j 
NEXT  i 

{a}  =  INVERSE ( [K] )  *  {f> 


Table  2.3:  MCI  Partial  Differential  Equation  Solution  Algorithm 


2.1.4  Conclusions 

The  Monte  Carlo  Integration  method  is  very  easy  to  implement  in  a  code.  However  it  suffers 
from  three  difficulties: 


1.  Because  values  are  obtained  from  averaging  the  neighbors,  derivatives  near  the  bound¬ 
aries  of  the  material  will  not  be  correct.  The  MCI  method  inherently  requires  a  sym¬ 
metric  distribution  of  nodes  around  each  particle  to  correctly  calculate  derivatives. 
When  nodes  are  vacant  on  a  given  side  of  the  particle,  the  derivatives  will  be  poorly  cal¬ 


culated.  There  are  various  ways  of  dealing  with  this  problem  (Bonet  and  Kulasegaram, 


2002),  however  the  complexity  is  more  severe  than  claimed. 


2.  Although  values  of  functions  can  be  obtained  by  multiplying  each  particle  by  a  “mass,” 
different  sized  particles  produce  values  which  are  not  correct.  The  MCI  method  works 
best  with  many  small  particles,  and  larger  particles  will  not  be  evaluated  correctly. 


3.  Since  the  method  is  based  off  of  averaging  values  from  neighbors,  there  is  an  inherent 
viscosity  associated  with  the  method.  We  could  reduce  the  averaging  of  values  by 
making  the  value  h  very  small  so  as  to  use  less  particles.  In  fact,  the  kernel  reproduces 
the  exact  value  when  h  — >  0.  However,  when  we  reduce  h  to  conserve  energy  in  the 
system,  we  lose  resolution.  When  we  lose  resolution,  we  lose  stability,  especially  in  the 
calculation  of  derivatives.  Hence,  there  is  a  contradiction. 
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2.2  Radial  Basis  Functions 


Instead  of  using  averages  from  neighbors,  it  is  possible  to  calculate  interpolants  by  first 
solving  a  system  of  equations  for  known  values.  By  doing  this,  we  can  ensure  behavior  of 
the  Dirac  delta  function  at  each  node.  Furthermore,  by  applying  differencing  operators  to 
the  system  of  equations,  we  can  also  calculate  derivatives  of  the  interpolants.  Problems 
associated  with  the  geometric  distribution  of  points  in  the  MCI  method  can  be  alleviated, 
however,  problems  of  computation  may  arise. 


2.2.1  Formulation 

We  say  that  our  continuous  function,  /,  is  a  linear  combination  of  a  basis  function,  0,  and 
an  unknown  weighting  value,  A  (Cheng  et  al.,  2003). 


n— 1 

f(x)  = 

3=0 


(2.13) 


To  perform  interpolation,  we  specify  that  f(x  =  xj  =  f(xt).  Using  this  prescription,  we 
solve  eqn(2.13)  to  find  the  weights,  A.  Because  all  values  of  the  weighting  functions  are 
nodally  based,  derivatives  may  pass  directly  to  the  basis  function. 


n—  1 

V/( x)  =  V  VA(r,)A( 

3=0 


Xn 


(2.14) 


The  laplacian  is  treated  in  the  same  way. 


Tl—1 


V2/($)  =  Y.  v>.WAfe,) 

3=0 


(2.15) 


Once  the  weights,  A,  are  known,  derivatives  can  be  found  simply  by  evaluating  eqn(2.14). 


2.2.2  Prescriptions  on  basis  functions 


Two  types  of  basis  functions  may  be  chosen.  Local  basis  functions  are  similar  to  those  used 
in  the  MCI  techniques.  Types  of  local  basis  functions  include  Gaussian  curves  or  inverse 
multiquadrics.  They  result  in  sparse  matrices,  but  have  extremely  poor  convergence  on  the 


boundaries.  Equation  (2.16)  shows  a  Gaussian  basis  function. 


4>  (r„)  =  e 


(2.16) 


Global  basis  functions  have  the  form  that  distant  particles  have  more  influence  than  neigh¬ 
boring  particles.  Types  of  global  basis  functions  include  thin  plate  splines  and  multiquadrics. 
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While  global  basis  functions  mean  that  particles  have  more  influence  from  distant  particles 
than  they  do  from  neighboring  particles.  Equation  (2.17)  shows  a  thin  plate  spline  function. 


</>  (r.j)  =  rj,  log  r 


(2.17) 


Multiquadrics  can  be  shown  to  have  exponential  convergence  (Cheng  et  ah,  2003),  however, 


the  matrix  which  results  may  be  poorly  conditioned  and  very  full.  Multiquadrics  have  been 
found  in  Lagrangian  meshless  methods  to  produce  the  best  calculations.  The  parameter 
n  involved  in  multiquadrics  affects  the  type  of  kernel.  If  n  is  small,  the  kernel  becomes  a 
compact,  inverse  multiquadric.  If  n  is  large,  the  kernel  becomes  very  global.  The  kernel  for 
multiquadrics  is  (when  n  =  3): 


</>  (rtJ)  =  ( h 2  +  r2  +  r2) 


n— 3/2 


(2.18) 


The  gradient  is: 


V i/l>  (r,,)  =  (2 n  -  3)  rk  (h2  +  r2  +  r2)"  5/2 


(2.19) 


The  laplacian  is: 

V2</>  (rtJ)  =  (2 n  -  3)  (2 h2  +  (2 n  -  3)  r2)  (h2  +  r2  +  r2)n  7/2  (2.20) 


2.2.3  Evaluation 


We  may  use  Radial  Basis  Functions  to  solve  a  partial  differential  equation  via  collocation 
(Fedoseyev  et  al.,  2002).  We  will  use  matrix  notation  using  Einstein  summation.  For  the 
field  equation  we  have: 


V23  =  (V2<Mrg)  A, 

(2.21) 

For  the  flux  boundaries,  we  solve: 

V/j  =  (V0„)  Xj 

(2.22) 

For  Dirichlet  boundaries,  we  solve: 

ft  Aj  Aj 

(2.23) 

Table [2T| shows  the  initialization  routine  for  the  RBF  method.  The  algorithm  is  identical 
to  the  initiation  routine  for  the  MCI  method  in  that  the  same  sort  of  operators  are  created. 

Table  [23]  shows  how  derivatives  may  be  calculated  using  RBFs.  What  is  required  first  is 
the  weights  A.  Afterwards,  the  differential  operators  are  simple  to  evaluate. 

Table |2.6| gives  an  algorithm  of  how  a  boundary  value  problem  can  be  solved  with  Radial 
Basis  Functions.  This  algorithm,  like  that  for  the  MCI  method,  involves  solving  an  equation 

[K\{X}  =  {f}. 
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FOR  i  =  0  to  (n-1) 

FOR  j  =  0  to  (n-1) 
CALCULATE  r[i,j] 
CALCULATE  phi[i,j] 
CALCULATE  dphix[i,j] 
CALCULATE  dphiy[i,j] 
CALCULATE  ddphi[i,j] 
NEXT  j 
NEXT  i 


eqn(2.2) 


eqn(2. 


eqni 

eqni 

eqni 


(2 


(2 


(2 


18) 

19) 

19) 

20) 


Table  2.4:  RBF  Initialization  Algorithm 


{lambda}  =  INVERSE ( [phi] )  *  {f} 

FOR  i  =  0  to  (n-1) 

FOR  j  =  0  to  (n-1) 

dfdx[i]  =  dphix[i,j]  *  lambda  [j] 
dfdy[i]  =  dphiy[i,j]  *  lambda  [j] 
NEXT  j 
NEXT  i 


eqn(2.13) 


eqn(214) 

eqn(2.14) 


Table  2.5:  RBF  Derivative  Calculation  Algorithm 


FOR  i  in  the  field  equation 
FOR  j  =  0  to  (n-1) 

K[i, j]  =  ddphi [i, j] 

NEXT  j 
NEXT  i 

FOR  i  having  Dirichlet  conditions 
FOR  j  =  0  to  (n-1) 

K[i,  j]  =  phi  [i ,  j] 

NEXT  j 
NEXT  i 

FOR  i  having  Neumann  conditions 
FOR  j  =  0  to  (n-1) 

K[i, j]  =  dphix[i, j]  *  nx[i]  +  dphiy[i,j]  *  ny[i] 
NEXT  j 
NEXT  i 

{lambda}  =  INVERSE([K])  *  {f} 


Table  2.6:  RBF  Partial  Differential  Equation  Solution  Algorithm 
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2.2.4  Conclusions 

Radial  Basis  Functions  suffer  difficulties  in  the  form: 

1.  A  local  basis  will  not  properly  calculate  derivatives  on  the  boundaries.  This  is  the  same 
difficulty  as  found  with  the  MCI  method.  Therefore,  a  global  basis  function  must  be 
utilized. 

2.  Global  basis  functions  suffer  from  massive  ill-conditioning.  Conditioning  numbers 
0(17)  are  common.  Furthermore,  when  particles  begin  to  move  and  two  particles 
move  very  close  to  each  other,  the  matrix  will  become  singular. 

3.  Global  basis  functions  do  not  have  a  good  physical  interpretation.  A  particle  should 
have  more  influence  from  its  neighbor  than  it  does  from  the  furthest  particle. 
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2.3  Multi-dimensional  Lagrange  Interpolation  Polyno¬ 
mials 

The  difficulty  with  the  MCI  method  is  that  while  we  desire  great  accuracy  with  h  — *  0,  we 
also  desire  stability  and  resolution  when  h  is  very  large.  Different  sized  particles  can  not  be 
used  as  this  will  create  problems  with  particle  “vacancies.”  Furthermore,  derivatives  on  the 
boundaries  are  poorly  calculated.  However,  what  is  desirable  is  that  the  matrices  are  very 
sparse  and  are  not  expensive  to  calculate. 

Radial  Basis  Functions  alleviate  the  accuracy  problems  with  MCI  methods.  However, 
we  get  further  problems  in  that  an  extra  matrix  must  be  inverted  for  every  calculation.  Not 
only  is  the  matrix  full,  but  it  usually  has  a  very  poor  conditioning  number.  Finally,  the 
physical  interpretation  of  radial  basis  functions  is  that  distant  particles  have  more  influence 
on  a  particle  than  neighbor  particles.  This  can  lead  a  fluid  simulation  to  become  unstable. 

Instead,  a  new  method  is  devised  using  a  series  of  Lagrange  polynomial  interpolations. 
The  purpose  of  the  MLIP  method  is  to  extend  the  resolution  to  interpolate  for  many  particles, 
while  still  retaining  the  property  of  a  Dirac  delta  function.  Furthermore,  the  matrix  to  be 
inverted  maintains  the  properties  of  being  sparse  and  diagonally  dominant. 

2.3.1  Formulation 

Suppose  that  instead  of  individual  “particles”  or  “nodes”  of  fluid,  we  discretize  our  domain 
into  actual  finite  volumes.  To  represent  the  specified  partition  of  fluid,  we  use  a  set  of 
quadrature  points  within  each  volume. 


A  diagram  of  integration 

•  The  red  particle  is  center  of  the  integra¬ 
tion  volume,  Q 

•  The  white  particles  are  neighbor  parti¬ 
cles 

•  The  length  h  is  the  interaction  radius 

•  The  blue  particles  are  the  positions  for 
Gaussian  quadrature 


Figure  2.2:  MLIP  Integration 

Our  purpose  is  to  relate  the  deformations  of  each  volume  to  all  others.  The  Lagrangian 
interpolation  formula  expresses  a  one- dimensional  function  g(x)  as  a  polynomial  ln(x)  of 
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degree  n  in  a  linear  equation  with  a  second  function  f(x)  (Hildebrand,  1987). 

n—  1 

g(x)  =  lQ(x)f(x o)  +  h(x)f(xi )  H - h  ln(x)f(xn)  =  lk(x)f(xk )  (2.24) 


fc=0 


The  lagrangian  coefficient  function  4(x)  is  a  solution  of  this  system  of  linear  equations  which 
also  has  the  property  of  lt{x3)  —  StJ.  The  coefficient  function  can  be  expressed  by  the  monic 
polynomial: 


n(x)  =  (x  —  Xo)(x  —  X\)  •  •  •  (x  —  xn) 

(2.25) 

The  solution  for  the  coefficient  function  becomes: 

/  \  k—l 

77(X)  TT  X  ~  Xk 

Ux)  =  7 - — t — r  =  - 

(X  ~  Xt)TT,{Xt)  fJJ  X l  ~  Xk 

(2.26) 

Changing  variables  of  g(x)  =  f(x),  we  obtain: 

n—  1 

zoo  =  + E(x ) 

3=0 

(2.27) 

The  associated  error  with  this  interpolation  is: 

w  \  ,  ^/n+1(0 
e(i)=i(i)(»+1)! 

(2.28) 

If  we  now  desire  to  calculate  the  one-dimensional  integral,  we  can 
Gaussian  quadrature  as  a  function  of  a  weight  IT 

express  it  through  a 

p  71—1 

/=  /  w(x)/(x)da:  «  W  Wjf(x3) 

j=0 

(2.29) 

The  weight  IT  can  be  evaluated  as  the  weight  from  a  quadrature  point  w 
interpolation: 

(x)  and  the  Lagrange 

71—1 

Wj  =  y ^wt(x)k(x) 

1=0 

(2.30) 

Written  in  index  notation  and  using  Einstein’s  summation  convention, 

f=  w(x)f(x)dx  «  w.jljJi 

(2.31) 

in 


In  one-dimension,  the  integration  occurs  within  the  limits  of  each  particle’s  size.  However, 
particles  “interact”  with  other  particles  outside  of  the  integration  domain  via  the  interpo¬ 
lating  polynomial  1. 
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Taking  derivatives  of  the  Lagrange  polynomial,  we  get: 


dlt(x) 

dx 


m—  1 

E 

m=  1 
m/z 


l 


X  Xm 


A: — 1 

n 

k= 0 

k^i 


X  ~  Xk 


X%  X]% 


(2.32) 


d\{x ) 
dx2 


n—  1 

E 

n=l 

n^z 


l 


x  —  xn 


m—  1 

E 

m=l 

m^i^n 


l 


X  Xm 


n— 1 

n 

fc=0 


X  ~  Xk 


xt  xk 


(2.33) 


2.3.2  Multi-dimensional  evaluation 

In  multiple  dimensions,  we  perform  the  interpolation  for  each  quadrature  point  in  the  par¬ 
ticle’s  volume.  Since  the  Lagrangian  polynomial  is  inherently  a  one-dimensional  function, 
dimensional  “cuts”  must  be  performed.  For  each  quadrature  point,  in  two-dimensions,  we 
divide  the  surrounding  region  into  two  cuts  ( x  and  y).  We  include  or  exclude  points  in  a 
cut  based  on  the  angle  to  the  point  (eqn(|2.34l)).  Each  cut  uses  one- dimensional  Lagrangian 
interpolation. 

0l3  =  arctan(Xj  —  £  )  (2.34) 


A  diagram  for  multi-dimensional  interpolation 

•  The  blue  point  is  the  quadrature  point 
to  perform  the  interpolation 

•  The  green  particles  are  neighboring  vol¬ 
umes  to  be  interpolated 

•  The  white  particles  are  excluded  from 
the  directional  interpolation 

•  The  gray  sectors  represent  the  spaces  to 
be  interpolated 

•  The  red  squares  are  the  calculated  “in¬ 
terpolation  points” 


Figure  2.3:  Dimensional-cut  for  the  MLIP  method 

The  interpolating  polynomial  for  each  cut  passes  through  every  point  we  give  to  it. 
Therefore,  if  we  have  a  large  number  of  points,  the  order  of  the  polynomial  will  become 
extreme,  and  the  interpolation  will  become  chaotic.  Therefore,  we  will  limit  the  number  of 
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points  used  in  the  interpolation  to  be  parabolic.  These  points  will  be  called  “interpolation 
points.” 

For  each  dimensional-cut,  we  can  divide  the  area  into  3  sectors  to  provide  3  interpolation 
points.  We  can  calculate  the  position  of  the  sectors  by  constructing  a  histogram  for  each  cut. 
Each  sector  from  the  histogram  contains  a  set  of  points  which  influence  each  center.  The 
location  of  each  “interpolation  point”  is  determined  by  a  weighted  average  of  the  particles 
in  each  sector. 

Table  |2.7|  shows  the  algorithm  to  create  a  histogram  for  a  dimensional  cut.  In  brief,  we 
create  a  MCI  smoothing  technique  with  a  very  large  smoothing  radius.  By  multiplying  the 
weighting  kernel  by  a  penalty  vector  initially  containing  values  of  unity,  we  can  find  the  most 
“dense”  sector  for  each  cut.  After  we  find  the  most  dense  area,  we  modify  the  penalty  vector 
to  be  negative  at  that  location.  We  repeat  the  process  with  the  modified  penalty  vector  to 
find  the  next  most  dense  sector. 


function  makebins(x,y) 

FOR  i  =  0  to  n 
FOR  j  =  0  to  n 

r[i,j]  =  radius (x,y) 
w[i,j]  =  exp(-abs(r[i, j] )) 
NEXT  j 

penalty(i)  =  1. 

NEXT  i 


eqn(2.2) 


FOR  bin  =  0  to  2 

location_of _max  =  max(w*penalty) 
points [bin]  =  x(location_of _max) 
penalty [location_of_max]  =  -100 
NEXT  bin 

RETURN (points) 


Table  2.7:  MLIP  Sectoring  Algorithm 


Now  that  we  have  each  cut  divided  into  sectors,  we  need  to  interpolate  for  each  “interpo¬ 
lation  point”  of  each  sector.  We  can  interpolate  for  each  point  by  either  MCI  or  RBF.  RBFs 
have  been  shown  to  provide  the  most  accurate  interpolations  within  each  sector  of  each  cut. 
By  using  RBFs  for  each  sector,  we  obtain  the  influence  for  each  point  by  inverting  a  very 
small  (and  not  ill-conditioned)  RBF  interpolation  matrix.  Table  2.8  shows  the  algorithm 


for  performing  a  local  sector  interpolation.  The  inputs  for  the  function  are  a  matrix  of  the 
distances,  r  for  each  particle  in  a  sector  to  each  other,  and  a  vector  of  the  distance  rcenter 
from  each  particle  to  the  “interpolation  point”  of  the  sector. 

The  procedure  for  this  method  is  best  performed  in  a  particle- by-particle  basis.  For 
each  particle,  we  call  an  initialization  algorithm  and  then  multiply  the  resulting  Lagrange 
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function  weight (r,r_center) 

FOR  i  =  0  to  n 
FOR  j  =  0  to  n 

weight [i ,  j]  =  phi(r [i, j] ) 
NEXT  j 

wbin[i]  =  phi (r_center [i] ) 
NEXT  i 


eqn(2.18) 


eqn(2.18) 


w  =  wbin  *  INVERSE (weight) 
RETURN (w) 


Table  2.8:  MLIP  Local  Sector  Interpolation  Algorithm 


polynomial  vectors  by  the  sector  weight.  This  algorithm  for  each  quadrature  point  i  is  shown 
in  table  12.91 


FOR  j  =  0  to  (n-1) 
CALCULATE  r[j] 
CALCULATE  theta [j] 
NEXT  j 


eqn(2.2) 

eqn(2.34) 


FOR  k  =  0  to  dim 
CREATE  points [k] 

FOR  p  =  0  to  2 
CALCULATE  w[j,p] 

NEXT  p 

1  +=  Lagrange  (points  [k] )  *  w[j,p] 
dl[k]  =  dLagrange  (points  [k] )  *  w[j,p] 
ddl  +=  ddLagrange (points [k] )  *  w[j,p] 
NEXT  k 


algorithm  (|2.7[) 
algorithm  (|2.8[) 


eqn 

eqn 


(2.26) 


(2.32) 


eqn(2.33) 


Table  2.9:  MLIP  Initialization  Algorithm 


Besides  the  initiation  routine,  we  can  solve  an  equation  or  calculate  derivatives  in  the 


same  manner  as  with  the  MCI  method  using  algorithms  shown  in  Tables  2.2  and  2.3 


2.3.3  Conclusions 

1.  The  MLIP  method  has  the  desirable  property  of  having  a  relatively  large  interaction 
radius  with  other  particles,  but  preserves  the  Dirac  delta  property. 

2.  The  matrices  created  using  the  MLIP  method  are  sparse  and  diagonally  dominant. 
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This  is  extremely  beneficial  for  machine  storage  and  the  expense  to  invert. 

3.  The  use  of  directional  cuts  and  choosing  interpolating  points  has  the  potential  of  pro¬ 
ducing  spurious  and  incorrect  results.  However,  by  using  quadrature  points  for  each 
particle,  the  effect  of  any  incorrect  results  can  be  seriously  reduced. 

4.  There  is  a  considerable  amount  of  overhead  in  performing  the  initialization  routine. 
However,  by  using  a  language  such  as  C++,  the  sparsity  of  the  method  can  be  utilized 
in  a  particle-by-particle  manner  through  the  creation  of  a  particle  class. 
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Chapter  3 
Results 


To  evaluate  meshless  Lagrangian  methods,  there  are  two  components  which  must  be  tested. 
The  first  component  to  be  evaluated  is  the  accuracy  of  the  computational  methods.  We 
will  evaluate  the  accuracy  of  the  Monte  Carlo  Integration,  Radial  Basis  Functions,  and 
Multi-dimensional  Lagrange  Interpolating  Polynomials.  For  the  same  type  of  function,  the 
computational  methods  were  tested  on  different  types  of  particle  arrangements. 

The  second  component  of  the  meshless  Lagrangian  methods  to  be  evaluated  is  the  accu¬ 
racy  of  the  fluid  dynamic  models.  The  equations  of  motion  were  applied  to  different  types  of 
simulations.  We  will  evaluate  the  properties  of  the  Moving  Particle  Semi-implicit  equations 
as  well  as  the  Lagrange  Implicit  Fraction  Step  equations. 


3.1  Methods  of  Computation 


There  were  three  different  particle  distributions  used  for  testing  methods  of  computation: 

•  Full  structured  grid 

•  Sparse  structured  grid 

•  Unstructured  grid 


The  simulation  involved  using  a  function  /  from  eqn(3.2)  and  creating  an  interpolation, 


gradients,  and  Laplacian  operators  for  each  method.  The  function  from  eqn(3.2)  is  set  up  to 


be  similar  to  the  pressure  from  a  breaking  dam  problem.  It  is  a  solution  of  a  boundary  value 
problem  (eqn(|3.1|))  with  a  domain  (f2):  x  =  (l,xmax)  and  y  =  (1,  ymax)-  The  free  surface  is 
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considered  to  be  on  x  =  xmax  and  y  =  ymax •  The  walls  are  on  x  =  1  and  y  =  1. 


1)  V2/  =  0  in  fl 

2)  /  =  0  on  Td,  the  free  surface 

o\  df  n  , 

3)  — —  =  0  on  x  =  1 

ox 

df 

4)  —  =  —64  on  y  =  1 
5?/ 


4x2 

i  max 

a  —  16 - 


V^l  ftxfnax  '^‘rnaxVmax 


^ maxV < 


max  Umax 


(3.1) 


The  solution  to  this  boundary  value  problem  is  our  function  /: 

/  =  1  (o  (x  -  l)2  -  a  (xmax)2)  (128  (y  -  1)  -  a  (y  -  if  -  128 ymax  +  a  ( ymaxf )  (3.2) 

with  the  parameter  a  given  as  a  function  of  the  domain: 


(3.3) 


Figures  [3TTH3.4|  show  the  function  /  as  well  as  the  analytical  derivatives. 

The  discrete  interpolation  was  conducted  for  all  three  methods  of  computation. 

•  Monte  Carlo  Integration 

•  Radial  Basis  Functions 

•  Multi-dimensional  Lagrange  Interpolating  Polynomials 

Error  estimates  for  an  analytical  function  g  and  an  interpolated  value  g  are  given  by: 

EIL'c1  \g(x  =  xl)-g(xl)\ 


error 


E:=ob(*  =  20l 


(3.4) 
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Figure  3.1: 


Function  /  of  eqn(3.2) 


Figure  3.2:  V2/  of  eqn(3.2) 


Figure  3.4: 


§£  of  eqn(3.2) 
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3.1.1  Full  structured  grid 

The  structured  grid  involves  144  particles.  It  has  a  range  from  1  to  12  in  both  x  and  y  axes. 
Figure  [3~5]  shows  the  arrangement  of  particles. 


Figure  3.5:  Full  structured  grid 


For  all  results  shown  below,  the  analytical  values  are  shown  by  a  colored  surface.  The 
discrete  values  shown  numerically  are  represented  by  blue  circles.  When  a  mesh  for  the 
interpolated  values  can  be  created,  these  values  are  linked  together  by  a  transparent  surface. 
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A)  Monte  Carlo  Integration 


Figure  3.6  shows  the  MCI  interpolation  on  a  structured  grid  as  obtained  by  eqn(|2.9[).  Notice 
that  there  is  a  small  amount  of  error  between  the  true  value  and  the  interpolated  value.  This 
is  due  to  an  averaging  between  nodes.  Figure  |3.7|  shows  the  Laplacian  interpolation  on  a 


structured  grid  from  eqn(2.11).  There  is  a  considerable  amount  of  error  at  the  boundaries. 
However,  the  error  in  the  center  of  the  domain  is  almost  negligible. 


error  =  0.0100867 


error  =  1 .1892 
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Figure  3.6:  MCI:  /  on  a  full  structured  grid 


Figure  3.7:  MCI:  V2/  on  a  full  structured  grid 


Gradients  are  obtained  from  algorithm (|2.2|).  Figure  3.8  shows  the  gradient  in  the  x 
direction  on  a  structured  grid.  There  is  an  extreme  amount  of  error  on  x  =  1.  Figure  |T9| 
shows  the  gradient  in  the  y  direction  on  a  structured  grid.  The  same  extreme  error  is  also 
present  on  y  =  1. 


error  =  0.977996  error  =  0.98091 8 


Figure  3.8:  MCI:  §|  on  a  full  structured  grid  Figure  3.9:  MCI:  on  a  full  structured  grid 

In  terms  of  the  derivatives  in  the  MCI  method,  boundary  points  suffer  from  accuracy 
because  of  a  vacancy  of  points  on  a  given  side.  This  error  is  increased  when  the  value  of  the 
gradient  is  very  large. 
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B)  Radial  Basis  Functions 


Figure  3.10  shows  the  RBF  interpolation  on  a  full  structured  grid  as  obtained  by  eqn(2.13). 
There  is  no  error  for  the  interpolation  since  the  values  are  solved  for  via  Gaussian  -  elimina¬ 
tion.  Figure  [3~TT|  shows  the  Laplacian  interpolation  on  a  structured  grid.  There  is  a  small 
amount  of  error  on  the  boundaries,  but  it  is  still  very  good. 


error  =  7.03272e-1 4  error  _  q.03431  59 


Figure  3.10:  RBF:  /  on  a  full  structured  grid 


Figure  3.11:  RBF:  V2/  on  a  full  structured  grid 


Figures  |3.12|  and  |3.13|  show  the  gradients 
error  on  these  grids. 


in  the  x  and  y  directions.  There  is  almost  no 
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C)  Multi-dimensional  Lagrange  Interpolating  Polynomials 


Figure  |3.14  shows  the  interpolation  from  the  MLIP  method.  The  matrix  generated  by  the 
function  on  this  structured  grid  is  the  identity  matrix.  Figure  3.15  shows  the  Laplacian. 
There  is  very  little  error  in  the  Laplacian. 


error  =  3.7001 2e-1 2 


error  =  0.0144234 


Figure  3.14:  MLIP:  /  on  a  full  structured  grid 


Figure  3.15:  MLIP:  V2/  on  a  full  structured 
grid 


Figures  |3.16|  and  |3.17|  show  the  gradients 
error  on  these  grids. 


in  the  x  and  y  directions.  There  is  almost  no 


error  =  0.00238278 


error  =  0.000661944 


ro 


Figure  3.17:  MLIP:  ^  on  a  full  structured  grid 


Figure  3.16:  MLIP:  ^  on  a  full  structured  grid 
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3.1.2  Structured  sparse  grid 


The  sparse  grid  is  generated  via  a  quadtree  method  (Wang  et  ah,  1999).  The  quadtree 
method  is  set  up  to  search  for  particles  near  boundaries.  Particles  near  boundaries  are 
subdivided  until  a  suitable  mesh  division  is  created. 

The  quadtree  method  has  very  nice  characteristics  of  placing  more  particles  where  needed. 
Furthermore,  where  a  large  number  of  particles  are  unnecessary,  these  areas  are  left  to  only 
a  few  particles  having  very  large  masses.  There  were  only  70  particles  used  in  the  sparse 
grid. 


Figure  3.18:  Sparse  structured  grid 
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A)  Monte  Carlo  Integration 


Figure  |3.19|  shows  the  MCI  interpolation  on  a  sparse  grid.  Because  of  the  large  separation 
between  points,  the  interpolation  actually  improves  on  this  grid  due  to  less  averaging.  Figure 
3.20|  shows  the  Laplacian  interpolation  on  a  sparse  grid.  In  calculating  derivatives,  the 
averaging  of  method  breaks  down  for  “large  particles.”  There  is  a  tremendous  amount  of 
error,  and  the  method  is  entirely  unusable  for  this  grid  arrangement. 


error  =  0.00526696 


error  =  1 .34851 


Figure  3.19:  MCI:  /  on  a  sparse  structured  grid 


Figure  3.20:  MCI:  V2/  on  a  sparse  structured 
grid 


Figures  |3.21|  and  |3.22|  show  the  gradients.  Gradient  operations  in  the  MCI  method 
perform  poorly  between  differently  sized  particles. 


error  =  4.17991 


error  =  2.50889 


Figure  3.21:  MCI:  on  a  sparse  structured  grid 


Figure  3.22:  MCI:  ^  on  a  sparse  structured  grid 
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B)  Radial  Basis  Functions 


Figure  [3.23|  shows  the  RBF  interpolation  on  a  sparse  grid.  Problems  of  “large  particles”  as¬ 
sociated  with  the  MCI  method  are  fully  alleviated.  However,  because  of  the  large  differences 
in  the  sizes  of  particles,  the  eigenvalues  of  the  matrix  become  more  separated,  leading  to  a 
larger  conditioning  number  in  the  matrix  to  be  inverted  from  algorithm(|2.6[).  The  condition¬ 
ing  number  to  solve  the  boundary  value  problem  is  3(10)9.  Figure  |3.24|  shows  the  Laplacian 
interpolation  on  a  sparse  grid.  The  error  on  the  boundaries  remains,  but  is  still  acceptable. 


error  =  2.35228e-14  error  =  0.121716 


Figure  3.23:  RBF:  /  on  a  sparse  structured  grid  Figure  3.24:  RBF:  /  on  a  sparse  structured  grid 


Figures  |3.25|  and  |3.26|  show  the  gradients 
error  on  these  grids. 


in  the  x  and  y  directions.  There  is  almost  no 


error  =  0.0241449 


error  =  0.0226777 


Figure  3.25:  RBF:  on  a  sparse  structured  Figure  3.26:  RBF:  on  a  sparse  structured 

grid  grid 
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C)  Multi-dimensional  Lagrange  Interpolating  Polynomials 


The  MLIP  method  works  well  on  a  sparse  grid  work  well  while  also  alleviating  the  problems 
of  the  MCI  method.  Furthermore,  the  eigenvalues  in  the  resulting  matrix  do  not  have  the 
same  problems  as  associated  with  Radial  Basis  Functions.  The  conditioning  number  to  solve 
the  boundary  value  problem  is  only  214.  Figure  [3728]  shows  the  Laplacian  interpolation  on  a 
sparse  grid.  There  is  even  less  error  than  with  the  Radial  Basis  Functions,  but  small  errors 
remain  on  the  boundaries. 


error  =  3.08831  e-1 2 


Figure  3.27:  MLIP:  /  on  a  sparse  structured  grid 


error  =  0.0418167 


Figure  3.28:  MLIP:  V2/  on  a  sparse  structured 
grid 


Figures  |3.29|  and  |3.30|  show  the  gradients  in  the  x  and  y  directions.  There  is  almost  no 


error  in  the  calculation  of  gradients. 
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Figure  3.29:  MLIP:  |  ona  sparse  structured 
grid 


Figure  3.30:  MLIP:  ^  on  a  sparse  structured 
grid 
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3.1.3  Unstructured  grid 

This  unstructured  grid  is  set  to  simulate  action  after  particles  have  been  moved  slightly 
relative  to  an  initially  structured  grid  as  in  figure  |3.5[  Some  particles  move  very  close  to 
other  particles,  while  other  particles  move  further  apart. 


10 


12 


Figure  3.31:  Unstructured  grid 
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A)  Monte  Carlo  Integration 


Figure  |3.32|  shows  the  MCI  interpolation  on  a  structured  grid.  The  error  present  in  the 
structured  grid  is  the  same  error  in  the  unstructured  grid.  However,  the  amount  of  error 
does  not  change  dramatically  after  particles  have  been  moved. 


error  =  0.0165428 


error  =  1 .23527 


Figure  3.32:  MCI:  /  on  an  unstructured  grid 


Figure  3.33:  MCI:  V2/  on  an  unstructured  grid 


Figures  |3.34|  and  |3.35|  show  the  gradients  in  the  x  and  y  directions.  The  error  on  the 
boundaries  is  still  present,  however,  the  error  in  the  middle  of  the  mesh  changes  only  slightly. 


error  =  2.08028 


error  =  2.00006 
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B)  Radial  Basis  Functions 


Figure  |3.36|  shows  the  RBF  interpolation  on  a  disordered  grid.  There  is  no  error  in  this 
interpolation.  However,  figure  |3.37|  shows  the  Laplacian  interpolation  on  a  disordered  grid. 
Because  certain  particles  are  too  close  to  one- another,  the  calculations  can  break  down  and 
produce  very  spurious  results.  When  applied  to  a  fluid  code,  this  can  cause  the  simulation 
to  break  down. 


error  =  8.46033e-14 


error  =  13544.2 


Figure  3.36:  RBF:  /  on  an  unstructured  grid 


Figure  3.37:  RBF:  V2/  on  an  unstructured  grid 


Figures  |3.38|  and  |3.39|  show  the  gradients  in  the  x  and  y  directions.  The  calculation  of 
gradients  is  still  very  good. 


error  =  0.00217386 
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Figure  3.38:  RBF:  ^  on  an  unstructured  grid 


Figure  3.39:  RBF:  on  an  unstructured  grid 
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C)  Multi-dimensional  Lagrange  Interpolating  Polynomials 


The  interpolation,  figure  [3T0j  is  shown  to  be  exact,  while  the  Laplacian,  figure  |3.41|  shows 
the  smallest  degree  of  error  from  any  method. 


error  =  1 .4291  e-1 2 


error  =  0.0301842 


Figure  3.40:  MLIP:  /  on  an  unstructured  grid 


Figure  3.41:  MLIP:  V2/  on  an  unstructured  grid 


The  gradients,  shown  in  figures  |3.42|  and  |3.43|  show  excellent  accuracy  for  a  disordered 
mesh. 


error  =  0.00468268  error  =  0.00244536 
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3.2  Fluid  Codes 


For  a  meshless  Lagrangian  method,  we  may  use  a  combination  of  any  equation  of  motion  and 
method  of  computation.  In  total,  there  are  nine  different  combinations  we  can  use,  shown 
in  table  13.11 


Equations  of  Motion 

Method  of  Computation 

Spherical  Particle  Hydrodynamics  (SPH) 
Moving  Particle  Semi-implicit  (MPS) 
Lagrange  Implicit  Fraction  Step  (LIFS) 

Monte  Carlo  Integration  (MCI) 

Radial  Basis  Functions  (RBF) 

Multi-d  Lagrange  Interpolating  Polynomials  (MLIP) 

Table  3.1:  Equations  of  Motion  and  Methods  of  Computations 


The  codes  chosen  for  testing  some  of  the  theories  included: 

A  MPS  and  Monte  Carlo  Integration  (MPS  -  MCI) 

B  Lagrange  Implicit  Fraction  Step  and  Radial  Basis  Functions  (LIFS  -  RBF) 

The  SPH  equations  were  not  investigated  because  a  large  amount  of  research  has  already 
been  conducted  on  the  subject.  Instead,  it  was  decided  to  investigate  implicit  codes  which 
did  not  use  an  equation  of  state  to  calculate  the  pressure.  Multi-dimensional  Lagrange 
Interpolating  Polynomials  also  were  not  used  because  an  efficient  C+- 1-  algorithm  has  not 
been  written  yet. 

The  MPS  -  MCI  code  was  written  in  C++.  The  code  made  use  of  the  sparsity  of  the  MCI 
method,  and  its  evaluation  was  very  quick.  The  LIFS  -  RBF  code  was  written  in  Matlab. 
Because  of  Matlab’s  robust  matrix  inversion,  it  was  utilized  since  the  matrices  involved  in 
Radial  Basis  Functions  are  very  full  and  ill-conditioned. 

Two  types  of  experiments  were  performed  using  both  methods.  Each  experiment  was 
used  to  evaluate  the  accuracy  and  robustness  of  the  codes.  To  simplify  the  problem  as  well 
as  to  investigate  the  accuracy  of  the  simulations,  viscosity  was  set  to  zero.  This  was  done  to 
evaluate  if  there  was  any  numerical  damping  involved  in  the  computations. 


41 


3.2.1  Dam  Break 


The  dam  break  involves  a  water  column  of  height  H  and  length  xo  in  which  one  side  of  the 
column  is  released  at  time  t  =  0+.  This  simulation  allowed  for  the  evaluation  of  boundary 
conditions  as  well  as  a  comparison  to  experimental  data. 
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A)  Moving  Particle  Semi-Implicit  and  Monte  Carlo  Integration 


The  propagation  of  the  water  can  be  compared  to  published  results  (Shao  and  Lo,  2003), 


shown  in  figure  3.45,  English  units  were  used  in  the  simulation  with  g  =  32.2 ft/s2  and 
p  =  1.99  slugs/ ft3. 


0  0.5  1  1.5  2 

Non-dimensional  time  T  =  t(g/H)1  / 2 


Figure  3.45:  Comparison  of  dam  breaks 


Figure  |3.46|  shows  the  propagation  of  the  water  as  it  decays  from  its  initial  position. 
Identification  of  the  free  surface  particles  utilized  inherent  defects  in  Monte  Carlo  Integration. 
Particle  vacancies  were  identified  by  a  lack  of  symmetry  in  the  derivative  matrices,  dwx  and 
dwy.  If  the  lack  of  symmetry  was  over  a  given  limit,  the  particle  would  be  tagged  as  being 
on  the  free  surface. 

Figure  |3.47|  shows  the  water  column  hitting  another  wall  on  the  opposite  side.  As  the 
water  descended  from  the  opposite  wall,  a  small  vortex  was  created  in  the  bottom  right  corner 
of  the  fluid.  The  presence  of  this  vorticity  indicated  that  there  was  numerical  damping  in 
the  Monte  Carlo  integration  technique. 
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Figure  3.47:  Water  column  hitting  a  wall  using  MPS  -  MCI 


45 


isii 


B)  Lagrange  Implicit  Fraction  Step  method  and  Radial  Basis  Functions 


The  water  column  collapse  shows  promise  for  the  LIFS  method  with  the  inclusion  of  correct 
boundary  conditions  on  the  fluid.  Figure  |3.48  shows  the  pressure  distribution  (as  contour 
lines)  as  well  as  the  velocity  (as  vectors)  when  t  =  0+. 

However,  the  use  of  Radial  Basis  Functions  proved  to  be  both  slow  and  unstable.  Figure 


3.49  shows  the  subsequent  decay  of  the  water  column.  As  a  particle  would  approach  another 


particle,  the  calculations  would  degenerate  because  of  a  singular  matrix  in  Radial  Basis 
Functions. 
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Figure  3.48:  Initial  pressure  distribution  in  the  water  column  using  LIFS  -  RBF 
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Figure  3.49:  Propagation  of  water  column  using  LIFS  -  RBF 
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Finally,  the  identification  of  free  surface  particles  was  only  done  at  the  beginning  of  the 
simulation.  Because  a  particle  was  not  re-tagged  if  a  free-surface  particle  fell  into  the  fluid, 
a  free-surface  particle  would  not  obey  wall  boundaries.  Figure  [3750]  shows  the  consequences 
of  not  re-identifying  boundary  particles. 


Figure  3.50: 


Degeneration  of  LIFS  -  RBF  solution  due  to  boundary  condition  mistracking 
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3.2.2  Sloshing  Wave 


The  sloshing  wave  experiment  had  an  initial  condition  of  particles  arranged  in  a  sinusoidal 
fashion.  Because  the  fluid  was  inviscid,  the  oscillations  of  the  fluid  should  continue  in  forever. 
Furthermore,  in  a  given  period,  the  particle  arrangement  should  always  regenerate  the  initial 
condition. 

The  specific  shape  of  the  initial  domain  is  given  by  figure  |3.51|  In  the  initial  position,  we 
can  calculate  the  maximum  hydrostatic  at  (x,y)  =  (0,0)  as: 


p  =  pg  (*o  +  <?) 


(3.5) 


Figure  3.51:  Diagram  of  Sloshing  Wave  Experiment 
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A)  Moving  Particle  Semi-Implicit  and  Monte  Carlo  Integration 


The  MPS  method  showed  that  there  is  a  large  amount  of  numerical  damping  in  the  sine- wave 
calculations.  The  simulation  also  required  an  unnecessary  amount  of  particles.  Because  a 
sparse  grid  could  not  be  utilized  with  the  Monte  Carlo  method,  particles  in  lower  portion 
of  the  domain  had  to  be  calculated,  but  had  little  consequence  on  the  free  surface.  These 
particles  remained  present  so  that  the  Monte  Carlo  method  would  not  break  down.  The  loss 


of  kinetic  energy  can  be  shown  in  figure  3.52 


time  (sec) 

Figure  3.52:  Loss  of  kinetic  energy  in  MPS-MCI  in  the  sloshing  wave  experiment 


Figures  |3.53|  and  |3.54|  show  the  first  period,  T  of 
tions,  small  vortices  were  seen  in  the  bottom  corners 


the  oscillations, 
of  the  domain. 


In  much  later  oscilla- 
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Figure  3.53:  MPS-MCI  Sloshing  solution  at  t/T  =  0  (top),  =  1/4  (bottom) 
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Figure  3.54:  MPS-MCI  Sloshing  solution  at  t/T  =  1/2  (top),  =  3/4  (bottom) 
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B)  Lagrange  Implicit  Fraction  Step  method  and  Radial  Basis  Functions 


The  sloshing  simulation  shows  that  the  algorithm  from  the  LIFS  method  is  viable,  however, 
the  use  of  Radial  Basis  Functions  in  large  deformation  problems  is  not.  The  pressure  and 


velocity  vectors  at  time  t  =  0+  are  shown  in  the  upper  figure  3.56  The  enforcement  of  flux 
boundary  conditions  in  the  partial  differential  equation  results  in  a  correct  calculation  of 
pressure. 


Pressure  at  (x,  y)  =  (0, 0)  can  be  compared  between  both  methods  in  figure  3.55.  Because 
flux  boundary  conditions  are  not  enforced,  the  MPS  equations  need  a  number  of  iterations  to 
correctly  calculate  the  pressure.  This  occurs  when  a  fluid  particle  “falls”  into  a  wall  particle 
to  such  a  point  that  the  field  equation  can  correctly  reproduce  the  required  boundary  force. 

The  propagation  of  the  sinusoid  is  seen  in  figures  |3.56| 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35  0.4 

Time  (sec) 


Figure  3.55:  Comparison  of  pressure  normalized  by  hydrostatic  pressure  eqn(3.5) 
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Figure  3.56:  Propagation  of  the  sinusoid  using  Radial  Basis  Functions 
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However,  after  large  deformations,  the  calculations  involved  with  RBFs  create  the  situa¬ 
tion  that  particles  in  the  held  equation  come  too  close  to  each  other.  The  problem  with  this 


is  that  the  RBF  calculations  become  unstable,  and  the  simulation  breaks  down.  Figure  [T57 
shows  the  creation  of  a  singularity  in  pressure  near  (x,y)  =  (60,40). 


Figure  3.57:  Degeneration  of  solution  due  to  singularity  in  Radial  Basis  Functions 
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3.3  Conclusions 


3.3.1  Methods  of  Computation 

The  Monte  Carlo  Integration  method  posed  two  sets  of  problems.  First,  the  calculation 
of  derivatives,  especially  on  the  boundaries  produced  very  poor  results.  Mathematically, 
derivatives  in  the  MCI  method  require  an  infinite  domain.  If  we  have  a  vacancy  on  any 
portion  of  the  domain,  the  calculation  of  derivatives  will  not  be  accurate. 

As  a  consequence,  flux  boundary  conditions  from  the  MCI  method  cannot  be  used  to 
solve  a  system  of  equations.  Therefore,  the  LIFS  method  is  not  usable  with  Monte  Carlo 
Integration. 

Secondly,  the  MCI  method  needs  a  constant  spacing  of  particles  to  obtain  accurate  deriva¬ 
tives.  This  prohibits  the  use  of  differently  sized  particles  in  a  simulation.  By  not  being  able 
to  use  differently  sized  particles,  the  application  of  the  MCI  method  for  any  real  naval  hy¬ 
drodynamics  or  full  three-dimensional  aerodynamics  problems  is  severely  limited. 

Radial  Basis  Functions  alleviate  many  of  the  geometric  problems  associated  with  the 
MCI  method.  However,  the  machine  costs  associated  with  inverting  a  full  and  ill-conditioned 
matrix  limits  the  amount  of  particles  which  could  be  used  in  calculations.  The  cost  to  invert 
an  RBF  matrix  can  never  be  less  than  0(n3).  This  cost  also  prohibits  the  application  of 
RBFs  for  any  simulation  with  a  large  number  of  particles. 

Furthermore,  if  a  grid  deformed  too  much,  the  matrix  could  become  singular,  and  the 
simulation  would  break  down.  It  is  recommended  that  Radial  Basis  Functions  be  used  in  an 
Eulerian  type  simulation  where  a  matrix  could  be  inverted  only  once  per  simulation. 

Multi-dimensional  Lagrange  Interpolating  Polynomials  were  created  to  preserve  the  nice 
machine  characteristics  of  the  MCI  method,  while  performing  more  accurate  interpolation. 
The  result  is  a  method  which  was  shown  to  provide  accurate  interpolation,  however,  it  is 
relatively  untested.  The  machine  expense  of  performing  local-interpolation  for  each  sector 
is  justified  by  the  very  low  conditioning  numbers  of  the  resulting  matrices  needed  to  solve  a 
boundary  value  problem. 

3.3.2  Fluid  Codes 

The  MPS  equations  combined  with  the  MCI  method  was  shown  to  provide  a  stable  and 
robust  technique.  However,  the  lack  of  a  flux  boundary  condition  for  the  differential  equation 
meant  that  the  problem  was  ill-posed.  Pressure  calculations  needed  multiple  time  steps  to 
allow  for  a  fluid  particle  to  “enter”  a  wall  particle.  Once  the  force  from  the  wall-particle 
sufficiently  repelled  the  fluid  particle  enough,  the  pressure  calculation  would  be  satisfied. 
Although  the  differential  equation  is  ill-posed,  the  MPS-MCI  code  still  works  because  the 
field  equation  enforced  on  the  wall  boundary  particles  reaches  an  equilibrium  after  enough 
time  steps. 

Furthermore,  in  the  inviscid  calculations,  small  vortices  would  occur  as  a  result  of  smooth¬ 
ing  from  the  MCI  technique.  Ultimately,  this  would  lead  to  a  decay  of  energy  in  the  system. 
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The  LIFS  equations  combined  with  RBFs  showed  that  the  inclusion  of  flux  boundary 
conditions  meant  that  every  time  step  was  capable  of  correctly  calculating  the  pressure-field. 
In  fact,  the  LIFS  equations  differed  from  MPS  only  in  the  initial  formulation  of  the  velocity 
components.  By  including  the  velocity  components  in  the  calculation  of  pressure,  a  well- 
posed  boundary  value  problem  was  formed.  Furthermore,  by  separating  the  influence  of 
viscosity  from  the  other  forces,  it  was  possible  to  form  a  more  efficient  inviscid  algorithm. 

3.3.3  Recommendations  and  Future  Work 

Boundary  conditions  still  pose  a  problem  for  meshless  Lagrangian  fluid  dynamics.  The 
identification  of  a  particle  being  either  on  the  free  surface  or  the  wall  needs  to  occur  for 
each  time  step.  If  not,  a  free  surface  particle  is  capable  of  “falling-through”  a  wall,  since  its 
pressure  is  always  set  to  zero.  The  identification  of  boundary  particles  on  a  sparse,  meshless 
grid  is  not  trivial.  One  possible  solution  could  be  a  quadtree  “remeshing”  of  the  domain 
after  a  certain  amount  of  deformation.  While  this  would  be  an  efficient  solution  which  is 
also  capable  of  preserving  boundaries,  this  violates  a  certain  trait  of  the  formulation  in  being 
meshless. 

A  second  problem  is  in  the  free-surface  particles.  Because  all  free-surface  particles  have 
zero  pressure,  the  pressure  gradient  perpendicular  to  the  normal  of  the  free  surface  is  always 
zero.  For  the  breaking  dam  experiment,  this  means  that  the  free-surface  particles  were  not 
capable  of  supporting  each  other  from  the  influence  of  gravity.  Instead,  they  will  always  fall 
into  each  other.  We  could  realize  that  the  pressure  for  a  free  surface  volume  is  not  actually 
zero  within  the  volume,  but  only  on  the  outer  surface  of  that  volume.  A  possible  solution 
could  include  creating  “phantom”  free  surface  particles  in  an  outward  normal  direction  one 
unit  length  away.  These  phantom  particles  would  have  zero  pressure,  allowing  the  “free- 
surface”  particles  to  have  a  pressure  able  to  support  one  another. 

A  third  problem  developed  where  a  particle  was  both  against  the  wall  and  on  the  free 
surface.  In  terms  of  solving  the  system  of  equations,  a  resolution  could  include  using  the 
flux  boundary  conditions  in  the  system,  but  while  also  applying  a  penalty  term  in  case  the 
pressure  was  not  zero.  This  problem  could  also  be  resolved  by  the  creation  of  the  “phantom 
particles”  for  the  free  surface. 
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Appendix  A 
Algorithms 


The  appendix  is  divided  as  follows:  Appendix  |A|  details  the  program  flow  of  the  MPS-MCI 
and  LIFS-RBF  codes.  Appendix  [B]  details  the  inputs  and  outputs  to  the  codes.  Appendix 
□  shows  the  codes  contained  on  the  CD. 
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A.l  MPS-MCI  Algorithm 
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A. 2  LIFS-RBF  Algorithm 
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Appendix  B 
Code  Explanation 


B.l  MPS-MCI  Algorithm 

The  MPS-MCI  code  takes  an  ASCII  input  file  for  arguments  to  the  domain  and  environ¬ 
mental  variables.  An  example  file  is  shown  as  dam-break. input: 

-32  .01  4 

1  0  2  1  1  10  20 

2  0  2  0  0  30  0 

2  0  2  0  1  0  30 

2  0  2  30  1  30  30 

On  the  first  line,  the  first  number  is  the  acceleration  due  to  gravity.  The  second  number  is 
the  value  of  St.  The  third  number  is  the  number  of  “domains.” 

The  subsequent  lines  are  inputs  to  the  “domains.”  The  first  number  qualifies  if  the 
domain  is  water  (1)  or  a  wall  (2).  The  second  number  is  a  “modification”  to  the  domain. 
The  modification  does  nothing  when  mod=0.  The  modification  changes  the  domain  into  a 
sloshing  experiment  when  mod=l.  The  third  number  is  the  density  of  each  domain. 

The  following  four  numbers  in  the  line  detail  the  geometry  of  each  domain.  The  numbers 
are  inputed  as  xO,  yO,  xl,  yl. 

The  output  from  the  program  are  two  files.  The  created  file,  out.avi ,  is  a  movie  of  the 
simulation.  The  second  created  file,  output,  is  an  ASCII  file  showing  the  maximum  velocity, 
maximum  pressure,  and  kinetic  energy  of  the  system  for  every  time  step. 

B.2  LIFS-RBF  Algorithm 

The  inputs  to  the  LIFS-RBF  code  are  written  directly  into  the  file  lifsrbf.m.  The  variables  of 
St,  gravitational  acceleration,  and  density  may  be  modified  in  line  11  of  the  file.  To  perform 
a  sloshing  experiment,  the  lines  17-25  should  be  uncommented.  To  perform  a  dam- breaking 
experiment,  the  lines  29-36  should  be  uncommented.  Output  is  sent  to  a  variable  mov.  This 
can  be  turned  into  an  animation  by  using  Matlab’s  built-in  movie  to  avi  function. 
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Appendix  C 
File  Listing 


The  CD  included  with  this  thesis  contains  the  MPS-MCI  code,  the  LIFS-RBF  code,  the 
methods  of  computation  algorithms,  and  this  document.  Presently,  this  CD  is  also  posted 
at:  www.redcoopers.net 

The  folder  “MPSMCF  contains  the  MPS-MCI  code: 

1.  main.cpp,  and  classes.  *  are  the  main  computational  functions  for  the  MPS-MCI  algo¬ 
rithm.  The  Numerical  Recipes  toolkit  is  used  for  mathematical  functions  and  arrays. 

This  toolkit  is  available  for  purchase  from  www  .nr.com.  The  GLUT  library  is  also  used 
to  display  the  simulations  to  the  screen. 

2.  avi*.  *,  opengl.  *,  and  utility.  *  are  hies  used  to  generate  movies  of  the  simulations. 

The  folder  “LIFSRBF”  contains  the  LIFS-RBF  code: 

1.  lifsrbf.m  is  the  main  code. 

2.  rbf-kernel.m  and  radius. m  are  functions  to  compute  the  RBF  operators. 

3.  sparseinput.m,  octree,  to,  and  uniquer.m  are  functions  to  create  a  sparse  mesh. 

4.  printf.  *  are  utilities  written  in  C  for  Matlab  for  printing  to  the  screen,  printf.mexmac 
must  be  recompiled  for  a  machine. 

The  folder  “MOC”  contains  the  methods  of  computation  testing  algorithms: 

1.  grad-2d.m  is  the  main  code. 

2.  radius,  m  calculates  the  distances  between  every  particle. 

3.  mcLkernel.m  is  the  code  to  generate  MCI  operators. 

4.  rbf-kernel.m  is  the  code  to  generate  RBF  operators. 

5.  lip -kernel,  m  is  the  main  function  to  generate  MLIP  operators.  It  uses  the  hies  liprbf-kernel.  to, 
makebins.m,  and  Lagrange*. m  as  utilities. 

6.  sparseinput.m ,  octree. to,  and  uniquer.m  are  functions  to  create  a  sparse  mesh. 
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